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Abstract 

An analytic representation of the short-range repulsion energy in ionic systems is described that 
allows for the fact that ions may change their size and shape depending on their environment. 
This function is extremely efficient to evaluate relative to previous methods of modeling the same 
physical effects. Using a well-defined parametrization procedure we have obtained parameter sets 
for this energy function that reproduce closely the density functional theory potential energy surface 
of bulk MgO. We show how excellent agreement can be obtained with experimental measurements 
of phonon frequencies and temperature and pressure dependences of the density by using this 
effective potential in conjunction with ab initio parametrization. 
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I. INTRODUCTION 

The problem of modelling the dynamics of ionic materials has a long history P, For 
many years the only available models were empirical force fields that were not very accurate 
or transferable between different environments. This was a problem both of the form of 
the model potentials used and the way in which these potentials were parametrized. Most 
frequently, simple pairwise effective potentials were used whose parameters were obtained 
using a combination of physical reasoning and empiricism. Such potentials could, at best, 
only be expected to produce qualitative or poorhj-quantitative results. 

The advent of ab initio molecular dynamics^, |4| (MD) brought about a dramatic im- 
provement in the accuracy with which the potential energy surface of the ions could be 
calculated, but this came with the price of an enormous increase in computational expense. 
In molecular dynamics simulations the precision with which thermodynamic properties can 
be calculated depends on the size of the system studied and the length of the simulation over 
which averages may be taken. For ab initio molecular dynamics one is generally confined to 
systems of around 100 atoms and simulation times of ~ 10 picoseconds and so the precision 
with which many properties may be calculated is poor. In addition, for highly viscous liq- 
uids (such as silica^, (|) or very harmonic crystals the timescales available within ab initio 
MD may not be sufficient to adequately equilibrate the systemp]]. Some ab initio-sized sys- 
tems may suffer from finite size effects. For all these reasons there are many applications 
(good examples being the determination of melting temperatures or finite-temperature elas- 
tic constants) that are very difficult to tackle with ab initio MD and others (such as thermal 
conductivity and viscosity) that cannot be addressed at all. 

Effective potentials are far less computationally expensive and so allow much longer simu- 
lation times and larger system sizes. It is therefore extremely desirable to find a compromise 
between the accuracy of ab initio MD and the computational speed of MD using effective 
potentials. 



It has been shown 



lOl lll | how the accuracy of an effective potential can be greatly 



improved by parametrizing using information from ab initio calculations. However in order 
for this to work a physically appropriate form for the potential must be used. Improving 
effective potentials therefore involves both using good parametrization procedures and ap- 
propriate forms for the potentials. The form that an effective potential takes should be 



flexible enough to describe the potential energy surface and specific enough to allow for 
efficient parametrization and application. 

In this paper we present a many-body force field for ionic systems that incorporates the 
effect of an ion's environment on its shape and size and the impact that such ionic distortions 
have on the short-range repulsive interactions between ions. These effects have been shown 
in the pa 8t to be i m po rt ant to the bo„di„ g of m a„ y simple ion.e SyS te mS Q UQ. The 
force field presented has much in common with previously proposed "compressible-ion" 3] 
and "aspherical-ion" 3] models but is superior from a computational point of view and 
avoids some problems associated with the way in which these models are implemented (see 
section ITTTj) . It can also be implemented within a general mathematical framework that 
is very amenable to changes for the purposes of experimentation and improvement. A 
particularly attractive feature of the proposed force-field is that it can describe distortions 
of arbitrary shape. Previous models have been confined to describing distortions of certain 
given symmetries. 

We apply the proposed model to crystalline and liquid magnesium oxide. MgO has been 



used in the past as a testing ground for models of ionic systems \14\. It is considered to 



be the simplest oxide and yet it is a system in which many-body interactions associated with 
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id ]. It is therefore a 



distortions of the large oxide anion are known to be important 
starting point for attempts to model the many body interactions in oxides and ionic systems 
in general. 

MgO is a system of geophysical importance. It is an important component of the earth's 



lower mantle and its stability up to high pressures 1 171] make it useful as a pressure calibration 
standard for high pressure and temperature experiments. Since the vast majority of the 
compounds present in the earth's mantle are oxides, MgO is a starting point for studies of 
the effects of pressure and temperature on mantle materials. 

For MgO the use of our model is found to result in a significant improvement (over 
pairwise representations of the short-range repulsion) in the ability to fit the forces, stresses, 
and energies extracted from density functional theory (DFT) [18] calculations. We use the 
model in conjunction with electrostatic interactions (involving both charges and induced 
dipoles) within a well-defined parametrization procedure. The resulting force field is shown 
to give a very good description of the structure and dynamics of MgO. 



II. PARAMETRIZING A FORCE FIELD USING AB INITIO DATA 



As mentioned above, it has been shown for a number of systems that a very high level of 
accuracy may be achieved by parametrizing an effective potential by fitting to DFT forces, 
stresses, and energies in selected atomic configurations However empiricism 

is also necessary unless one uses a functional form that is physically appropriate in the 
sense that it describes the electronic effects that are most important for the property one 
is studying j]]]. In particular in reference Q we have shown that, for silica, the inadequacy 
of simple pair potentials means that improving agreement with the forces and stresses from 
ab initio simulations, in general, disimproves the ability of empirically derived parameter 
sets to describe crystal structures. By inclusion of more many-body effects (in this case the 
polarizability of the oxygen anion) the ability of the force-field to reproduce ab initio data is 
improved and the empiricism becomes unnecessary for describing structure very accurately. 

Here we make the assumption that ab initio calculations are highly accurate and always 
superior to calculations using effective potentials. This assumption is clearly formally un- 
justifiable but it is useful for our purposes. We aspire to ab initio accuracy and leave the 
improvement of the ab initio calculations themselves as a separate problem. We assume 
that if one achieves a perfect fit to ab initio data one gets an extremely accurate effective 
potential. However, as one moves away from this limit, the relationship between the fit to 
the ab initio data and the quality of the effective potential clearly weakens. If the fit is very 
poor then (as with the description of silica using pair potentials p|) even large improvements 
to it may not improve the ability of the potential to describe physical properties. 

We define the fit as 7 = (AF, AS, AE), where AF,AS and AE are dimensionless quan- 
tities corresponding to the average |12| percentage differences between forces, stress compo- 
nents, and energy differences between different configurations calculated ab initio and with 
the effective potential, i.e. 



AF = 100 x 



AS = 100 x 
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where Ff is the a th cartesian component of the forces on ion J, S al3 denote stress components, 
U k is the potential energy of configuration k, {77} is the set of parameters that, along with 
the functional form, characterise the effective potential. B is a pressure that may be taken 
to be the bulk modulus of the material. The superscripts 'ep' and 'ai' indicate whether a 
quantity has been calculated with the effective potential or ab initio. In order for 7 to be 
meaningful it is necessary that it be converged with respect to the number of configurations 
n c that have been used to test the effective potential. Typical values of n c are of the order 
of 10. 

Since our parametrization scheme does not discriminate between different contributions 
to the DFT forces, in order to be sure a priori that a potential works well one must have 
a very high quality fit to the ab initio data. This is because different properties depend to 
varying degrees on different contributions to the forces on the ions. It is difficult to know 
how a given value of, for example, AF, manifests itself in the property that one is interested 
in studying. This is particularly true when there are a number of important contributions 
to the forces as these different contributions may be fit by an effective potential to varying 
degrees. If one has a potential that fits the ab initio forces perfectly then, clearly, each 
contribution is fit perfectly and this uncertainty is eliminated. For example, for silica the 
structure was extremely well described by using a good description of electrostatic effects ?!]. 
However the forces differed from the ab initio ones by ~ 16%. This suggests that the gross 
features of the potential energy surface are well described but that more local features and 
finer details, that are important for dynamics, may not be. We have observed that the 
diffusion of liquid silica at 3500 K as described by the model proposed in reference jl seems 
to be too fast compared to extrapolations of lower temperature experimental data 20(. We 
also observe the phase transition between a— and f3— cristobalite at too low a temperature. 
Both of these observations are consistent with an underestimation of local energy barriers 
due to a functional form that oversimplifies the short-ranged interactions between ions. 

For every given functional form there are minimum values of AF,AS, and AE in param- 
eter space. Functional forms that oversimplify the description of the electronic effects that 
are relevant for a given physical property may not be capable of achieving a close enough 



fit to the ab initio data to ensure that improvement of this fit improves the ability of the 
parameter set to describe physical properties. If this is the case then it is still conceivable 
that one may improve the agreement with experiment on many physical properties in an 
empirical way We are of the opinion that an inability of a functional form to fit ab initio 
data indicates an intrinsic inadequacy of this form and that potentials created in this way 
should not be relied upon as being predictive. Effective potentials represent electronic ef- 
fects in a phenomenological way. A potential that cannot produce very realistic values of the 
forces, stresses and energy differences must lack an important ingredient in its functional 
form and even if it produces good values for certain physical properties (such as those to 
which it was fit), for other properties it may show qualitatively different behaviour to the 
real system. A good potential form should fit the ab initio data to a high degree and signif- 
icantly improving this fit should improve the property that one is interested in simulating. 
If this is the case then one can be confident that the ability to describe the property is due 
to a good microscopic description of the interatomic interactions. 

Ultimately, our aim is to produce potentials that may be relied upon quantitatively to 
predict, not only the structures and the thermodynamic properties of bulk ionic systems, 
but also their dynamical properties. Although much work has been done in this direction 
(see, for example, references 15 andQ and references therein) this is still a very ambitious 
goal. Finding a potential that describes dynamics well is a particularly difficult task ^| . 

However, the fact that DFT calculations can provide an essentially limitless amount 
of information that can be used in the parametrization process means that we are not 
constrained to using potential forms either with a small number of parameters or with 
parameters that have an obvious physical interpretation. 

We use the same iterative parametrization scheme that was used in reference Q. This is a 
slightly adapted form of the schemes used in references lol and llfll and it allows us to make the 
very specific and non-trivial statement about a potential that for any atomic configuration 
created with the potential under specified thermodynamic conditions, the forces are, on 
average, within pf% of those calculated ab initio 2jJ, the stress components within p s %, and 
the energy-differences between configurations within p e %. 

The parametrization scheme involves minimizing the function T({r]}) = ufAF + ljsAS + 
ujeAE with respect to the set of parameters {rj}. In order to be sure that the minimization 
procedure is meaningful the number of ab initio configurations used in the fit ( n c ) is 



required to be reasonably large. In general its value depends on the system studied, the 
potential form used and the number of atoms in the unit cell. It was found that for MgO a 
value of n c = 10 was required in order to achieve convergence in the fit 7. In each step of 
the iteration at least a further five configurations were retained during the fitting procedure 
in order to test that the final functional form fit these configurations as well as it did those 
that were used in the minimization of rfjr?}). 

n 

Minimization was performed using a combination of simulated annealing 22 and Powell 
. A basin in the surface defined by T({r]}) in 77— space was initially found 
using simulated annealing and, once found, further minimization was performed using the 
method of Powell. Minimization of r({r]})using realistic force fields, and particularly simu- 
lated annealing, is a very computationally expensive process. However simulated annealing 
is very useful for two reasons : 1. It is very stable; Powell minimization can break down if 
numerical errors (such as overflow errors) occur due to unphysical values of the parameters; 
2. In principle it can always bring one to the global minimum. In practice however this 
depends on how much computer time one is willing to allocate it. These properties of the 
simulated annealing method make it particularly useful when fitting a potential for the first 
time. One does not need to start with reasonable or physical values of the parameters in 
order for it to converge and this means that one may parametrize exotic potentials for which 
the parameters have no obvious physical interpretation. 

The freedom that one is afforded using the combination of ab initio data and simulated 
annealing is crucial. It simply would not be possible to parametrize a force field such as the 
distortable-ion model introduced in section IIVI without either one of these assets. If one 
depends on the use of experimental data, only force-fields containing a very small number 
of parameters may be parametrized if the fit is to be statistically significant. The amount 
of information that can be extracted from ab initio calculations is very large by comparison 
and so it allows much more complicated force-fields to be parametrized. If parametrization 
of a potential depends on good initial guesses for the parameters (as is the case with most 
optimization algorithms), each parameter must have an obvious physical interpretation and 
this can limit ones freedom when constructing a potential form. This problem can be solved 
by using simulated annealing in the parametrization process. 

The ab initio calculations were performed with DFT within the local density 



approximation 



2J, |2j| using the planewave pseudopotential method 26, 27]. We have used 



soft norm-conserving pseudopotentials[28] that are identical to those that have previously 



been used success 
and temperatures 



ully to calculate its vibrational properties for a large range of pressures 
2fl| . We require our simulations to produce good quality forces and Karki 
et al. have performed much more rigorous tests of these pseudopotentials than it would 
have been feasible for us to perform. We perform our calculations on unit cells containing 
64 atoms under periodic boundary conditions and the Brillouin zone is sampled using only 
the T— point. A kinetic energy cutoff of 120 Ryd. was used in the plane wave expansion of 
the wavefunctions in order to converge the stress very well. 

III. MANY-BODY MODELS FOR IONIC SYSTEMS 

In this section a system for including many-body interactions in ionic systems will be 
introduced. We confine our attention to compounds in which variations in the degree of 
ionicity (either locally or globally) may be neglected. For such systems many body con- 
tributions to the interatomic interactions arise from the response of the size and shape of 
the ions to their environment. The anions carry multipole moments and, particularly the 
lowest order of these, the dipole moments, have an important impact on phonon spectra. 
The change in size and shape of the anions also has an impact on the short ranged Pauli 
exclusion repulsion between ions. 
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14. 



15 



Although these effects are all interdependent, it has recently been shown 
how an improved description of ionic dynamics may be achieved by independently incorpo- 
rating them in the interatomic potential. Some of the ingredients of this force field have 
recently been parametrized using forces and stresses from density functional theory calcula- 
tions and this resulted in an extremely good description of the potential energy surface in 



the crystal 



111- 



The non-electrostatic repulsion between ions at short distances has in the past usually 
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been modelled as an exponential function of the interionic distance[l|, |3lJ. The use of this 
exponential form rests on the assumptions that ions are spherical and of fixed size and 
that the repulsion between them due to the Pauli exclusion principle is proportional to the 
overlap between the ions whose charge distribution tails off exponentially. Although this 
may be an adequate approximation in crystals of high symmetry at low temperatures and a 
given pressure, at higher temperatures or at a different pressure or when a change of phase 



occurs, it is likely that anions will readjust their size and shape to fill the available space. 

In order to cope with this, and in an attempt to improve the ability of ionic models 
to reproduce experimental equations of state and the relative energetics of different crystal 
structures, Wilson et al. have developed a compressible ion (CI) model Previously, 



ionic breathing effects had generally been incorporated within the shell model^^(for a 
recent application see reference l32l) . 

Wilson et al. wrote the potential energy due to short-range repulsive interactions between 
anion and cation as 

V^_({Rj}- {5aj}) = V scli ({5aj}) + Vo^Rj}; {Sen}) (1) 

where V^ e if is the sum of the changes in the internal energies of the ions and V ov is the total 
potential energy due to the repulsive overlap interaction between the ions, and 5ai is the 
change in the radius of ion / from its average value oj. From electronic structure calculations 
of the perfect cubic crystal it was found that V^ e if could be written as 

V scli ({Rj}; {6aj}) = £ D lC osh((36aj) (2) 

Ie- 

and the standard exponential form was adopted for the interaction energy : 

V or ({R I };{Sa I })= ]T B Ije - a ^-^ +s ^ (3) 

Je-,Je+ 

At each timestep, during a simulation, the {Saj} were required to take values that minimized 
V sv . 

Although this has been successful in reproducing some low temperature properties of 
MgO and CaO, such as the crystal energies as a function of volume, the prediction of phonon 
frequencies with this model was found to be quite poor^J. This is because the distortion 



of the anions are in general not spherically symmetric if the crystal is disordered [l^. Il5j|. To 
account for aspherical distortions of the anion, Rowley et al. have extended the previous 
compressible-ion model by introducing eight further degrees of freedom to model distortions 
of dipolar and quadrupolar symmetry. Self-energy functions associated with these degrees 
of freedom were postulated. 

With the aspherical-ion (Al) model and including dispersion effects and polarization 
effects, Rowley et al. managed to find parameters that gave good phonon dispersion curves 



for MgO. Very recently Aguado et al. have used an almost identical model, but with many 
of the parameters found from density functional theory calculations, to produce phonon 
dispersion curves and thermal expansion curves of a very high quality[ll]. Since thermal 
expansion depends on the second derivatives of the potential energy with respect to the ionic 
positions, this is a good test of this model's representation of the potential energy surface 
of the crystal. 

This method has the drawback that distortions of the anions are restricted to those of 
dipolar and quadrupolar symmetry. Although distortions of different symmetry could be 
incorporated, this would involve a large increase in the number of degrees of freedom and a 
corresponding decrease in the efficiency of the method. It would be desirable to include the 
effect of distortions of arbitrary shape in a way that is practicable. 

In order to implement the above model in a computationally efficient way, an extended ' 
grangian approach has generally been used that is analogous to the Car-Parrinello method 
for electronic structure theory. The degrees of freedom that describe ionic distortions ({Sai} 
in the simple compressible-ion case) are the variables with which a "fictitious" dynamics was 
associated. 

There are some problems associated with the use of a Car-Parrinello approach to mod- 
elling ionic distortions however: 



a- 
3, 



Recent work 33] has shown that, particularly for disordered systems, this method 
results in systematic errors in dynamics relative to methods in which the fictitious 
variables take their minimum energy values. 

The use of a small timestep is intrinsic to the extended lagrangian approach. For 



example, in the work of Aguado et al. ll| a step of only 0.05 femtoseconds was used. 
This is more than an order of magnitude smaller than the time step that can typically 
be used for non-Car-Parrinello approaches. 

The extended lagrangian approach is considerably more efficient during dynamics when 
variables are evolved from previous time steps than it is when minimizing the energy 
with respect to the extended variables for the first time. Since the parametrization 
process involves a large number of such minimizations, and is computationally very 
expensive, the extended lagrangian approach becomes much less efficient. 



IV. A DISTORTABLE ION MODEL 



Our goal is to develop a general mathematical framework within which the many-body 
effects of ionic compression and aspherical ionic distortion can be modelled, but that avoids 
the problems associated with the previously proposed models. In particular we would like 
to find a computationally efficient way of modelling these effects that avoids the use of an 
extended lagrangian and the introduction of fictitious degrees of freedom. 

There are many open questions regarding the interactions between simple (electronically 
isolated) ions and so the framework should be general enough to allow for continual im- 
provement of the specific form of the potential. For example, the form of the self-energy 
in equation El and the self energies associated with aspherical distortions can certainly be 
improved. It may also be worth investigating different forms of the overlap energy. 

We will be primarily concerned with the anion-cation interaction. Of much lesser concern, 
initially at least, is the anion-anion interaction energy that has been found to provide only 
~ 3% of the energetics of the perfect crystal jl^l Il5| . Although, the same cannot be said 
with any degree of certainty of more disordered phases, or systems of different stoichiometry 
such as Si02, it is nonetheless the most obvious place to start when constructing a potential. 

We assume that a distortable ion ( such as 2 ~ ) has its shape and size "influenced" 
by all sufficiently close neighbouring ions. Much as in the scheme of Wilson, Madden and 
coworkers [lZ ll^. an ion is described as a nucleus surrounded by a single membrane (rep- 
resenting the electrons) the radius of which is allowed to vary with the two polar angles 
(although in their case, the radius only varied in certain symmetric ways). The influence an 
ion J exerts on ion I can loosely be thought of as a restraining force on the ion's tendency 
to expand and this restraint has a dependence on the polar angles (9, &) in the spherical 
coordinate system centered on /. We also assume that the influence exerted at coordinates 
(9, 6) is zero if the angle between the outward unit vector at those coordinates £(9, 6) and 
the vector Rjj = Rj — R/ is greater than 90°. 

We write the total influence on I at (9, &) due to all the other ions as 

pf\o, <j>) = J2 MRuW, 0) • xjje(£(0, 0) • XJ/ ) (4) 

where Xjj = R ji/ Rji and 

Q(£(9.d) ■jc.tt) = 1 if £(9.6) -x (f > and if 1(9. 6) ■ x., r < (5) 



Functional forms for fu(R) will be discussed in Section IV. A. The subscripts IJ are to 
indicate that a different function is used for each distinct pair of ionic species. 
Apart from a multiplicative constant, the spherical average of pf > (9,4>) is 

P?> = (6) 

We write the angular dependent radius , <Ji(9, 0) of an ion as 

^M) = 4V o) )+-iVM>,0)) (7) 

In other words, the radius at (6>, 0) is written as a sum of an average value due to the 
influence of all the ions and deviations from that average. Functional forms for af^ and af^ 
will be discussed in the next section. The distance between the membranes of ions / and J 
along their line of centers is 

Lu = Ru - oi{d JU - oj^/j, <p u ) (8) 

where 6jj and 0j/ are defined such that 1(6 jj, = xj/. We will use the notation 

P?j = P?(0ji,<t>Ji) (9) 
4) = 4\p {O) ,P?(0ji,<f>Ji)) (10) 
07j = <Ji{0ji,(t>Ji) (11) 

We now define the contribution to the total energy of the system from the short-range 
repulsive interactions as a sum of pairwise interactions between membranes. 

u SR = £ uTAL^giARu) (12) 
i,j>i 

where gu(R) takes the value 1 for R < R a , for R > Rf, and decays smoothly from 1 to 
between R a and Rb. This allows us to truncate the interaction at intermediate distances. 
The expression for the forces on the ions is given in Appendix A. 

A. Applying the model 

In order to apply this model we clearly need to find suitable expressions for the functions 
fr r,cri \ and ai 1 } . We begin bv making the assumption that the most imoortant interaction 



is the anion-cation interaction although this will be extended at a later stage to include the 
anion-anion interaction in a limited way. For the moment we are concerned with systems 
with two species such as MgO and we assume that the cation is small and rigid. For MgO 
it is likely that this is a very good assumption, given its degree of ionicity. 

In order to draw a correspondence with the compressible ion model of Wilson et al. ^| 
(see section IIII|) we write the total energy of the system due to short-range repulsion as 



It- It-,Je+ 

+ e~ a — RlJ + B ++e~ a++RlJ (13) 

I,Je-,J>I I,Je+,J>I 

The values of the anion radii at any time should be such that this repulsive energy is 
minimized. In other words 

QySK 



daf 



, V/ (14) 



da) 



Je 



To simplify the notation we write B' = a^ + B + _ and ((a) 



(o)\ _ dV & 



do j 



C(af ) e - Q - + ^ 0> = - J2 B'e- a - +RlJ (16) 



As was noted previously in section HVl there has been some discussion about the form of 
the self-energy of compressible ions. Although in the paper by Wilson et ai.fl^ the form 
used was that of a hyperbolic cosine of the amount of compression Sa, applications of the 
shell-model generally use a harmonic expression for this energy and in a very recent 

paper [34 1 Marks et al. have argued that for the oxide ion one should treat the 2p 6 shell and 
the s 2 shells separately with harmonic and exponential compression energies respectively. In 
references Q and 



3J justification of the forms of V^ e if used were based on quantum chemical 



calculations of the cold crystal under a large range pressures. The self-energy of ions in 
disordered systems may be quite different to that in the perfect crystal and anyway, for 
interatomic distances that one should expect to find in the crystal at zero pressure and 
temperatures up to the melting point, the calculated self-energy is still close to a linear 
regime. Furthermore we do not see any compelling physical reasoning behind any of the 
forms used. For these reasons we consider the form of this enerev to be an open Question. 



As a preliminary test of our model we have chosen an exponential form for V^ e if as this 
simplifies the mathematics considerably. C( <J i°' > ) will then also have an exponential form so 
we may write 

A ie -^ e- a - + ^ = -J2 B'e~ a -+ R " (17) 

By merging constant terms to simplify the notation, this equation can be rewritten in 
the form 

af\pf) = C 1 + C 2 Hpf) (18) 

where we say that 

pf = Y,C*e c *» (19) 
j 

and Ci,C2..etc are constants. By analogy with equation |H1 we can say that 

fu = C 3 e~ c * R " (20) 

One is not confined to such simple forms for the self-energy but for many forms one cannot 
write equation El in terms of af^ and one is forced to find crf^ by an iterative procedure. 
This only has a very slight impact on the efficiency of calculating the potential. Another 
form that we have tried, and for which this procedure is used is 

, (o) + / , v» 

where 61,62-. etc are constants. This form was chosen according to the (admittedly, highly 
simplistic) physical reasoning that the internal factors that determine an ion's radius are 
the electrostatic energy which varies like the inverse of a distance and the kinetic energy of 
the electrons which varies like the inverse of a distance squared. 

The above analysis has shown that the distortable-ion model presented is mathematically 
equivalent to the compressible-ion model of Wilson et al. if af] = in the limit that the 
fictitious mass of the extended-lagrangian approach goes to zero. 

It is not possible to map our approach onto the aspherical-ion model. However, we take 
a different approach to aspherical distortions. Given the functions fjj and starting; 
point we may postulate a form for erf]. We assume that the local distortion of an ion's 
membrane scales with the local "density" in the same way as the average radius scales 

with the spherical average of the density p^ . We therefore write 

(i) 

4)=C^{%) (22) 



Since we will be parametrizing this force by fitting to ab initio data, the minimization 
routine has the freedom either to make the constant C5 very small or zero if this is not 
a reasonable functional form, or if aspherical distortions are really energetically equivalent 
to spherical ones it can make C2 = C5 in which case the distortions of purely spherical 
symmetry disappear and 



Although all the above derivation has assumed that this potential is only to be used 
for modelling cation-anion interactions, the generality and freedom afforded us by our 
parametrization procedure means that we lose nothing by trying to apply it to the anion- 
anion interaction also. We have done this by fitting parameters for the anion- anion inter- 
action and we have found that it does improve the ability of the model to fit the ab initio 
forces. A more sensible, but also more expensive way of tackling the anion-anion interaction 
would be to introduce a self-consistent procedure to minimize the angular dependent radii 
simultaneously. 

We also note that, as has been pointed out by Marks et al, different electronic shells have 
different compression characteristics. This could be modelled within the present scheme by 
having two or more such distortable ion potentials acting in parallel. 

In our application of this model we have generally used timesteps of between 1 fs and 1.5 
fs and the model has been found to be extremely efficient. As an example: in a simulation of 
512 atoms of crystalline MgO at 3300 K (using potential F which is discussed in section ITVC|) 
we have used a time step of 1 fs. In this simulation, the time required for each calculation 
of the distortable-ion contribution to energy, forces, and stress was 0.57 seconds on a single 
300Mhz SGI origin MIPS R12000 processor. On the other hand, the contribution of all 
the electrostatic terms was 3.0 seconds. The total energy in this simulation drifted by 
approximately 2 K per picosecond during this simulation. This energy drift can be almost 
eliminated by more fully converging the polarization at each time step. 

A point to note regarding equation 0]is that there is a discontinuous change in O when 
£(6, 0) • xj/ = 0. Although equation 0] itself is not discontinuous, its derivative with respect 
to the ionic positions is. This leads to discontinuous forces and consequently to a drift in the 
total energy. In practice we have not found this to be a problem as this drift in energy is very 
small compared to the drift that is due to the incomplete convergence of the polarization. 
However, if necessarv this problem mav be eliminated bv replacing 6 with a function that 




(23) 



varies smoothly from 1 to as £ ■ Xjj approaches zero. 



B. Testing the model 

The model that we have proposed is considerably more complex than a pair potential 
and so in testing the model we first would like to verify that it improves upon simpler force 
fields. We confine ourselves to testing the usefulness of three different ingredients of the 
interionic potential: 

1. A pairwise short-range interaction potential of the form 

where Bjj, ajj, Cjj, Ejj, Njj are all parameters to be optimized. 

2. A polarizable-ion potential including short-range polarization, as discussed in refer- 
ence Q- Only the oxygen ion is considered polarizable. 

3. A distortable-ion potential, as discussed in sections IIVI and II V Al The interaction 
energy between ions I and J is given by 

Uf^Ljj) = A Ije - aijLlJ + B IJ e~ PljLlJ (25) 

and the functions fu,cr\ , and crjj are given the same forms as in equations HH1 HE] 
and l2*2l respectively, i.e. 

fu = C?)e- C % R » (26) 
af\ P f)=CrH P f) (27) 

4}=C?]ln(% (28) 

where C\ from equation 1221 has been merged into the pre-exponential factors Ajj and 
Bij. The parameters to be optimized are Au, ajj, Bu, flu, Cj 1 ], Cf), Cj ,CjJ. 
The values R a = 8.5 a.u and Rb = 10 a.u. were used in the decay function gu. 

In addition to the ingredients mentioned, each force field also included the point-charge 
electrostatic potential with the charge on an ion as a parameter. 
We have parametrized five different force fields, as follows: 



A. 



A pair-potential : short-range pair potential, parametrized in the crystal at ambient 
conditions. 



B. A polarizable potential : short-range pair potential, with polarizable anions, 
parametrized in the crystal at ambient conditions. 

C. A distortable-ion potential : distortable-ion potential, parametrized in the crystal at 
ambient conditions. 

D. The full model : distortable-ion potential, with polarizable anions, parametrized in 
the crystal at ambient conditions. 

E. The full model : distortable-ion potential, with polarizable anions, parametrized in 
the liquid at 3000 K. 

In order to avoid the computational expense of performing a full self-consistent 
parametrization procedure for each of these potentials we have only used this procedure to 
parametrize three potentials using the full model (distortable-ion model with point charges 
and dipole polarization). The three potentials were optimised at zero pressure for (i) the 
liquid at 3000 K (ii) the crystal at 2000 K and (hi) the crystal at 300 K respectively. Each 
of these three potentials was used to create trajectories at the conditions for which it was 
optimised from which 'snapshot' atomic configurations were extracted and used in ab initio 
calculations. We will show that the full model is the best at reproducing the ab initio data 
and so these configurations were considered to be as realistic as we have the ability to create. 

Each of the potentials A to E was then parametrized using ab initio data from 10 of these 
configurations (at the relevant conditions). Each potential thus created was then tested 
on its ability to fit 10 new configurations (i.e. that were not used in the parametrization 
process) and also 10 configurations at each of the other two sets of conditions. For example, 
potential A was parametrized at 300 K and then its ability to fit the ab initio data in the 
crystal at 2000 K and the liquid at 3000 K was also tested. In all cases the error in the 
stress was evaluated relative to a pressure of B = 140 GPa. 

The results are summarized in table HI We cannot guarantee that we have found the global 
minimum in each case during optimization as simulated annealing had to be done at a rather 



rapid quench rate. The simulated annealing was followed by Powell minimization 23]. In 



each case, the total minimization time was the same (10 davs on a single processor") and 



therefore more economical force-fields are likely to be better minimized than less economical 
ones. A number of things can clearly be seen from table HI First of all, not surprisingly, 
the distortable-ion model on its own is quite bad. This is probably because of the shortness 
of the range of its interactions. Ions further away from each other than 10 a.u. interact 
only via the coulomb force between their charges. At 300 K, the full model is clearly better 
than all other forms. It also transfers very well up to higher temperatures and to the liquid. 
The pair-potential, although working quite well for the crystal, does not transfer well to 
the liquid. The polarizable model yields results that are intermediate in quality between 
the pair-potential and the full model. The results are a clear illustration of the fact that 
by adding more physics into the form of an effective potential one can create force-fields 
with, not only an improved ability to fit the ab initio data, but also a much improved 
transferability between different phases and conditions. 

The poor fit of potential E to the energy differences in the crystal at ambient conditions 
is because the energy differences in the liquid and high temperature solid are much greater 
than those at lower temperatures. The absolute value of the error in the energy differences 
is the same at low temperature and high temperature but AE is the error relative to the 
root-mean-squared value, which for the crystal is very small. 



1. Phonon Frequencies 

Having established that our inclusion of many-body effects has improved the potential 
form with respect to the pair potential, at least according to the criterion that we have 
adopted, we now look at its ability to model the vibrational spectrum of MgO. We note 
once again that the DFT scheme to which the potential was fit gives a very good description 
of phonon frequencies at ambient conditions [29]. 

We find the phonon frequencies for our potential at a number of k-points from the po- 
sitions of the peaks in the spectra of the spatial Fourier components of longitudinal and 
transverse charge and mass current correlation functions)^. 

We performed an MD simulation on a system of 512 atoms using the full-model, optimised 
in the crystal at 300 K (potential D). The current correlation functions were calculated on 
a time domain of length 2.9 ps that was averaged over a simulation of length 20 ps. The 
phonon dispersions that we get are shown in figure ^ We get an extremely close fit to both 



the experimental and the self-consistent DFT data. The chief discrepancies are in the optical 
modes which are systematically underestimated. The longitunical optical mode in particular 
is underestimated near the zone center. Although we do not calculate the mode frequencies 
at T = (0,0,0), as this would require an infinitely large simulation cell with the method 
that we are using, it looks as though the splitting between the longitudinal optical (LO) 
and the transverse optical (TO) phonons is slightly underestimated. In our parametrization 
procedure we have used a small cell to perform the ab initio calculations and so the long-range 
interactions that are important for dispersion near the T— point are not included. Our hope is 
that by modelling correctly the electrostatics at shorter range, we get a potential that, when 
used in a larger simulation cell, can accurately model the long range electrostatic interactions. 
This is not guaranteed however and is likely to work only if we include all relevant screening 
mechanisms in our functional form. The incorrect LO-TO phonon splitting suggests that 
our description of the electrostatics is incomplete since it arises from the long-range electric 
field induced by the LO phonon. This is not surprising since dipole polarization is only one 
of many screening mechanisms that are present in the real system. It may be that charge 
transfer between ions is important. However, a comparison with the results of reference [V. 
is suggestive of it being due to the fact that we haven't included the affects of higher-order 
multipoles. Density functional perturbation theory calculations [29J] of the Born effective 
charges yield values of Z* B = ±1.94. However, the charge on the oxygen ion in this potential 
(and all other potentials that we have fit) is ~ 1.5 - considerably less than this. Under 
the assumption that, within our model, short-range interactions and electrostatics describe 
completely separate aspects of the potential-energy surface (we do not know the extent to 
which this is true) the minimization routine fits the charge and the polarizability so as best to 
approximate the electrostatics of the crystal. The lack of higher order multipoles means that 
it must choose a compromise between purely dipole screening, in which the polarizability 
a and the charge q take their "true" values, and uniform screening in which the charge is 
simply reduced by a factor equal to the dielectric constant and the polarizability is zero. In 



reference 



11 



they use formal ionic charges and include both quadrupoles and dipoles and they 
get better agreement with experiment. Nevertheless, the description of the electrostatics that 
we have is significantly better than any other effective potential that we are aware of and 
quadrupoles would add considerably to the computational expense of the model. 



It is also worth noting that, although in reference II II the potential used seems to give a 



slightly better description of the optical phonon frequencies, they report a fit to the DFT 
forces of AF < 10%. Potential D is in significantly better agreement with our DFT data 
than this. Direct comparison of these two fits is not entirely justified, however. Aguado et 
al. fit to ab initio calculations of three different crystalline phases and therefore to three 
different coordination environments. Our fit of potential D is restricted to the phase and 
conditions at which we calculate the phonon frequencies. Although our fit is converged with 
respect to the number of atomic configurations tested, the small size of the fitting cell and 
the low temperature and high symmetry of these fitting conditions mean that only a very 
small region of phase space is visited. It may be that our functional form can numerically 
reproduce the potential energy surface in this small region of phase space in a way that is 
not necessarily physical and that therefore does not extend very well to describe the longer 
wavelength phonons that are present when we move to the larger simulation cell used for 
the calculation of phonon frequencies. 

On the other hand, if one is to construct a potential for use in a given region of phase 
space it may not always be advisable to include configurations from regions of phase space 
that are too far away from this domain of application. There is a danger that changes in 
electronic structure that cannot be represented phenomenologically by the functional form 
of the potential may occur. For example, a change in coordination of the ions may be 
accompanied by a change in the degree of ionicity. Since our model keeps the charges on 
each ion fixed, by fitting to different coordination environments we may find charges that 
are the best compromise between the optimal values for each different environment, and 
that are therefore ideal for none of them. For many materials, the liquid structure is not too 
dissimilar to the solid at the same pressure and temperature and may therefore be a good 
compromise as a fitting environment. By fitting to the liquid we can visit a large volume of 
phase space that is nevertheless centered on or near the crystalline phase of interest. 

We have used the above argument simply as a transparent way of discussing issues to 
consider when parametrizing a potential. The success of the potential created by Aguado et 
al. suggests that, for MgO, changes of ionicity do not present a problem in practice. 



C. Altering the model 



The simple exponential form that we have used for the self energy was chosen for its 
simplicity. We would like to test other forms of this energy and although we do not go in 
detail into this problem here, we do parametrize a potential using the form of equation PHI 
This has an appealing physical form and it requires one to self-consistently find the values 
of the radii af^ and af] . We observed for the previously created potentials that changing 
the radial cutoffs for the decay function gu of the distortable ion model, R a and Rf>, did not 
significantly change the fit to the ab initio. For larger distances the forces involved are either 
too small or cancel one another out. For this reason, and in order to improve efficiency we 
have used the slightly smaller values of R a = 7.0 a.u. and = 8.0 a.u. 

Using the full model we perform a full self consistent parametrization in the liquid at 3000 
K (potential F) and the solid at 3000 K (potential G). We have consistently found during 
our parametrizations that no significant improvement to the fit is obtained by including 
dispersion terms. It is likely that this is, at least partly, related to the limitations of the local 
density approximation as there is no reason to believe that it can describe such highly non- 
local dynamical fluctuations. Nevertheless, in the parametrization of potential F, 1/R 6 and 
1/R S terms were used in conjunction with Tang-Toennies dispersion damping functions [39!] . 
The fact that the parametrization procedure made some of these interactions repulsive (see 
table IHI)) is a clear indication that they should not necessarily be interpreted as representing 
the effect of electron dispersion. 

We have then tested the ability of potentials F and G to fit ab initio data at these 
conditions. For this testing we use 20 solid configurations and 20 liquid configurations. The 
results are summmarized in table |H] The fit is extremely good and even better than that 
obtained with the previous form of the self-energy. However, a problem has been encountered 
with the potential that was parametrized on the solid at 3000 K (potential G). We found 
that at higher temperatures and in the liquid, the iterative procedure by which the radii 
(o~u) were found sometimes failed to converge. This is the explanation for the poor fit of 
potential G to the liquid data. This highlights the importance of caution when applying a 
force field to conditions different from those in which it was parametrized. During extensive 
simulations of the liquid the distortable ion model for the potential parametrized on the 
liquid never failed to converge. 



The phonon dispersions for potential G are shown in figure 121 One should not expect 
results that are as good as those for a potential that is parametrized at ambient conditions, 
and so the results are extremely good. There is very good agreement with both experiment 
and the DFPT results of Karki et at. As before, the worst agreement is for the long- 
wavelength LO phonons, and once again this is probably due to our incomplete description 
of electronic screening. It may also be that the very high symmetry of the relatively cold 
crystal makes polarization energetically unfavourable, and so the polarizability appropriate 
for a hot crystal is larger. A too-large polarizability should manifest itself in the phonon 
curves as a lowering of the energy of the long-wavelength LO phonon modes due to improved 
screening of the macroscopic electric field. Nevertheless, in general the results seem even 
better than those of figure ^ and the ability of both potentials to reproduce ab initio energy 
differences is very satisfying and suggests that the form of the distortable-ion self-energy 
used may be better than a simple exponential. 



1. Density 

Figure El shows the equation of state of crystalline MgO at 300 K for potentials D and G. 
Although both are in extremely good agreement with experiment at low pressures, potential 
D is better at much higher pressures. There are a number of possible reasons for this. First of 
all, although we should not expect our DFT calculations to produce exactly the same equa- 



tion of state as that of Karki et a/|29| due to differences in the details of the calculations [35], 
their equation of state does become increasingly inaccurate at high pressures. It may be 
that the improved fit of potential G to the ab initio data results in a disimprovement in the 
equation of state due to an inadequacy of the ab initio data. Another possibility is that the 
reduced values of R a and Rb result in a reduction of transferability to high pressures due to 
the gradual introduction of more shells of neighbours in the distortable-ion calculation. 

Figure 13] shows the density as a function of temperature for potentials F and G (Potential 
D gave very similar results to potential G for a system size of 512 atoms). The experimentally 
observed thermal expansion is clearly very well reproduced by our potential. What is striking 
is that, particularly at low temperatures, finite size effects are very small. In reference 
finite size effects are much bigger. It may be that our model benefits from fitting the 
electrostatic interactions to the DFT data so that screening of interactions is more effective. 



It may also be that our potential over-screens the electrostatic interactions, i.e. that the 
dipole polarization is too large. This would be consistent with our underestimation of the 
LO mode frequencies approaching the T-point. 

Our potentials overall give an extremely good description of the density as a function of 
pressure and temperature. 

V. ANALYSING THE MODEL 

As discussed in section IIV Al we do not impose the distortable-ion model on the system. 
We have parametrized the force-field using simulated annealing that was begun at a high 
temperature. This means that, although we have supplied a functional form that is capa- 
ble of including distortable-ion behaviour, the minimization routine is free to do with this 
form whatever is best for reproducing ab initio forces. The options that are open to the 
minimization routine are 

• to disable all variable-radius functionality, and therefore to model the interionic forces 
with a double exponential of the interionic distance Rij. It would be optimal to do 
this if the way in which we model distortions is completely unphysical. 

• to enable only the compressible-ion part of the model, i.e. that which is analogous 
to the model of Wilson et a/Q], thereby allowing only spherically symmetric anion 
distortions. It would be optimal for it to do this if the way in which we model aspherical 
distortions is unphysical but our description of spherical distortions is reasonable. 

• to enable only the asymmetric part of the model and to disable purely spherically- 
symmetric distortions. This is optimal if our reasoning that aspherical distortions 
are energetically equivalent to spherical ones is true and the form of the model is 
reasonable. 

• to partially enable either or both types of distortions as the best compromise between 
rigid-ion behaviour, breathing-ion behaviour and distortable-ion behaviour if all three 
of the models fail to varying degrees and in different ways to reproduce the ab initio 
potential energy surface. 



The parametrization process is therefore itself a test of the distortable-ion model. We now 
look at what, precisely this parametrization process has done by examining the radius of 
an oxygen ion in the direction of a neighbouring magnesium ion for one of our potentials 
(potential F that is discussed in section HV C|) . The test is performed in the crystal at 3000K. 
The local radii of the anions consist of an arbitrary constant, that may be merged into the 
constant coefficient of the exponential force between ions, and the true variations of the 
radii due to changing environment. We look at the quantities ajj — WJJ , af] — crj-j, and 
o~f^ — o~f^ for anion I and cation J where aTj,af) and af^ are averages over a long trajectory. 
These quantities are therefore the non-constant parts of the different contributions to the 
radius of ion / in the direction of J (Recall that 07 j = + af] where includes only 
spherically-symmetric distortions and ajj includes aspherical distortions ). 

The results are shown in figure El and the variation in the value of Ljj, as defined by equa- 
tion |H1 along the same trajectory is shown for comparison. As can be seen, the local radius 
is dominated by the effect of the aspherical part of the distortable-ion model. The spherical 
part makes a significantly smaller contribution. This clearly vindicates our extension of the 
compressible-ion model to include aspherical distortions. 

The variation in the radius is very small compared to the variation in Ljj and so we look 
at what contribution this makes to the forces between the ions. Looking at the forces in a 
pairwise way is not entirely justified given the many-body nature of the potential, however 
it seems natural to look at the quantities 

dUffjLu) _ dUffiRu-an-oJi) 

Q?J = 100 x -J*U — - dR » (29) 

and 

dUff jLu) _ dUff jRu-WTJ-^JT) 

«S = 100 " M " «^<t,„" "' (30) 

dLu 

where Fj'f 1 ^' is the root mean-squared value (averaged over time) of the total force on 
anion / (i.e. from all atoms and from both electrostatic and non-electrostatic contributions) 
projected onto the line joining the centers of I and J. These quantities are plotted in figure^! 

is a way of looking at the impact of instantaneous variations of the membrane radii on 
the total force on the ion. is a way of looking at the impact of instantaneous variations 
of the membrane radii on just the short-range part of the force between ions I and J. If the 
radius of the ion is constant, then = = 0. It is difficult to know how one should 



best compare forces, or judge the impact of individual contributions to the forces. However, 
inspection of these two quantities strongly suggests that, with the parameters of the model 
chosen by the minimization routine, the variation of the anion's radius has a significant 
impact on dynamics. 

The above discussion shows that the minimization routine finds it optimal to allow fully 
aspherical distortions of the anions that impact significantly on the interatomic forces. 



VI. DISCUSSION 



In this paper we have presented a many-body interatomic potential for ionic systems that 
attempts to model the effect on interionic interactions of ions breathing and distorting as 
their environments change. The potential has been used in conjunction with a rigorous ab 
initio parametrization scheme and applied to the case of bulk magnesium oxide. The only 
empiricism involved in our construction of these potentials has been in the choice of the 
form of the potential. We have justified this form by its ability to fit the density functional 
theory potential energy surface and to reproduce experimental data. 

We have clearly shown that the form of the potential presented has a significantly im- 
proved ability (with respect to pairwise interactions) to fit the ionic potential energy surface 
of the hot crystal and the liquid calculated within density functional theory. It also has an 
improved transferability between phases and different physical conditions. 

Once parametrized carefully, the potential has been shown to produce excellent phonon 
dispersion curves, equations of state and thermal expansion. 

The computational expense of this potential is very small compared to other ingredients 
of realistic potentials such as dipole polarization. It is also li kely to be much less expensive 
than previously proposed extended lagrangian force- fields)^, 3], particularly since it allows 
the use of much larger timesteps. 

For all these reasons, the potential proposed is a valuable addition to effective force fields 
that aim towards a quantitative description of real ionic materials. 

The mathematical form of the potential is highly amenable to improvement and research 
into the optimal forms of its constituent functions may be expected to improve its ability to 
fit ab initio data. 



Appendix A 



Here we derive the expression for the forces on the ions from the definition of energy 
given in section HV1 

Using equations |H1 and the ath force component on ion K may be written as 

p«- V „ du "(¥*(x a \ daiJ dajf 

E UTj ^xUSik - Sj K ) (31) 
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where A IJK = 0(xjj • X/#). The notation Y^l(i) nas been introduced to indicate that 
the summation is over all ions L that are neighbours of I. This is necessary for practical 
implementation due to the truncation of interactions and to avoid summing over all the 
particles. 

Expanding equation EH we get 
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In the derivation we have made the further assumptions that fu = fji, Iff J = Ujf and 

9u = gji- 
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TABLES 



TABLE I: The fit to the ab initio data for the different potential forms 





300K crystal 


2000K crystal 


3000K Liquid 




AF 


AS 


AF 


AF 


AS 


AF 


AF 


AS 


AF 


A 


9.3 


5.0 


25.5 


13.7 


3.3 


15.5 


25.1 


4.8 


52.4 


B 


6.9 


5.2 


23.8 


9.0 


6.2 


17.8 


17.5 


5.6 


23.6 


C 


10.4 


39.1 


5.9 


13.6 


51.7 


164.8 


32.2 


58.2 


69.2 


D 


3.4 


0.6 


3.0 


6.8 


0.3 


9.8 


17.1 


0.3 


10.5 


E 


12.8 


0.1 


59.0 


10.2 


0.1 


18.9 


9.6 


0.0 


17.7 



TABLE II: The fit to the LDA ab initio data for the liquid (F) and solid (G) potentials used in 
the calculation of the melting slope. 





3000K Crystal 


3000K Liquid 




AF 


AS 


AF 


AF 


AS 


AF 


F 


9.6 


0.1 


10.8 


10.4 


0.2 


10.2 


G 


6.2 


0.3 


10.6 


44.0 


2.3 


54.0 



TABLE III: Parameters for potential F (atomic units) 



Parameter 


Mg-Mg 


Mg-0 


0-0 


b po \ a 


- 


1.55713 


4.01338 


Cpol 


- 


-1.28035 


31.93748 


C(l)6 




7.63066 x 10 1 


1.58612 x 10 3 


C(2)6 




2.01709 


2.42329 


(7(4) c 




-4.09954 x 10~ 2 


2.248341 x 10~ 2 


A d 


4.10638 x 10 12 


1.82404 x 10 2 


-2.78524 x 10 4 


a d 


7.29702 


2.22211 


2.98764 


B d 


-3.48800 x 10 13 


-2.00799 x 10 3 


4.00143 x 10 3 


13 d 


7.80852 


3.71852 


2.44883 



Parameter 


Mg 





Q 


1.44831 


-1.44831 


a 




14.25305 


C (3) 




-1.20698 x 10~ 2 






-1.44105 x 10 7 


g 2 / 




-5.14776 






1.11436 x 10 8 






7.20899 



Parameter 


Mg-Mg 


Mg-0 


0-0 






-1.01505 x 10 3 


3.24685 x 10 4 






9.16150 x 10 3 , -6.70301 x 10 6 


-5.93805 x 10 6 


b?t 9 




2.51113 x 10" 1 


-1.60446 x 10" 1 






3.78280 x 10- 1 , 9.07552 x 10~ 2 


-7.45819 x 10~ 2 







"Defined as in equation 4 of reference 

b See equation |2U 
c See equation 

d See equation 
e See equation 1771 
'See equation 1211 
^Defined as in Ref. 
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TABLE IV: Parameters for potential G (atomic units) 



Parameter 


Mg-Mg 


Mg-0 


0-0 


b P o\ 


- 


1.65744 


4.01338 


Cpol 


- 


-1.35136 


31.93748 


C(i) 




1.00281 


1.72179 x 10" 5 


C (2) 




0.65944 


0.19425 


C (4) 




-3.18047 x 10~ 2 


-2.93661 x 10~ 2 




1.39895 x 10 11 


2.16223 x 10 2 


-5.03851 x 10 4 


a 


6.72904 


2.23873 


3.09196 


B 


-2.50649 x 10 13 


-4.21223 x 10 2 


4.29645 x 10 3 


(3 


7.98372 


3.16438 


2.43995 



Parameter 


Mg 





Q 


1.48077 


-1.48077 


a 




10.40993 


C(3) 




-1.33660 x 10~ 2 


ei 




-5.59871 x 10 s 


£2 




7.61483 x lO" 1 


£3 




4.89669 x 10 8 


C4 




4.59914 x 10" 1 



Figures 



Fig. ^ The phonon dispersions of MgO as calculated with the full polarizable and 
distortable-ion model parametrized in the crystal under ambient conditions (Potential 
D) compared with experiment j3] and with the density functional perturbation theory 
results of Karki et al. 

Fig. El The phonon dispersions of MgO as calculated with the full polarizable and 
distortable-ion model parametrized in the crystal at 3000K (Potential G) compared 



with experiment |38l] and with the density functional perturbation theory results of 
Karki et a/.Q 



Fig. |3] The pressure of MgO as a function of density (at 300 K) compared to experiment 



MD simulations used the 



and density functional perturbation theory calculations 
full model potential parametrized at ambient conditions and simulation cells containing 
512 atoms. 

Fig. 13] The density of MgO as a function of temperature (at zero pressure) compared to 



experiment 



41 



42j for potentials G and F. Three different simulation cell sizes are used 



to check for finite size effects. 



Fig. El a) Ljj — Ljj as a function of time, where Ljj pa 4.66 a.u. is the average over the 
trajectory of Ljj, the inter-membrane distance (see section Hvjl. b) crij — oTJ, af] — afj , 
and <jj°^ — cxf^ . I and J are neighbouring anion and cation respectively. 

Fig. ElQ/j and QfJ (see equationsl2^1andl3l]|)as a function of time along the same trajectory 
shown in figure El 
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